Finite dissipation and intermittency in magnetohydrodynamics 
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We present an analysis of data stemming from numerical simulations of decaying magnetohydro- 
dynamic (MHD) turbulence up to grid resolution of 1536 3 points and up to Taylor Reynolds number 
of ~ 1200. The initial conditions are such that the initial velocity and magnetic fields are helical and 
in equipartition, while their correlation is negligible. Analyzing the data at the peak of dissipation, 
we show that the dissipation in MHD seems to asymptote to a constant as the Reynolds number 
increases, thereby strengthening the possibility of fast reconnection events in the solar environment 
for very large Reynolds numbers. Furthermore, intermittency of MHD flows, as determined by the 
spectrum of anomalous exponents of structure functions of the velocity and the magnetic field, is 
stronger than for fluids, confirming earlier results; however, we also find that there is a measurable 
difference between the exponents of the velocity and those of the magnetic field, as observed recently 
in the solar wind. Finally, we discuss the spectral scaling laws that arise in this flow. 
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As observations of astrophysical flows become more 
detailed, both spatially and temporally, the need for a 
deeper understanding of turbulent flows grows. In many 
such flows the fluid is coupled to a magnetic field the dy- 
namics of which can be understood in the magnetohydro- 
dynamic (MHD) approximation granted the analysis is 
confined to the large scales and under the hypothesis that 
the velocities are substantially smaller than the speed of 
light so that the displacement current in Maxwell's equa- 
tions can be neglected. This latter condition is easily 
fulfilled, characteristic bulk velocities in the solar wind 
being typically between 400 and 800 km s _1 , and tur- 
bulent velocities in the solar convection zone being ~ 1 
km s _1 . The MHD approach breaks down at small scales 
where kinetic plasma effects become important and one 
needs to include other terms in a generalized Ohm's law, 
such as ambipolar diffusion in weakly ionized plasmas as 
encountered in the interstellar medium, the Hall current 
for highly ionized media such as the solar wind, or an 
anisotropic pressure tensor. In such cases, the nonlincar- 
ities of the dynamical equations become more numerous 
and complex, parameter space is expanded and the re- 
sulting problem is quite challenging. For that reason, 
MHD is still a valid approach, albeit a simplified one, to 
tackle questions concerning the fate of a turbulent fluid 
coupled to magnetic fields. 

Laboratory experiments have classically been one 
venue to understand the physics of such fluids, for ex- 
ample in the context of reconnection [l[. Using liquid 
metals in the laboratory is a challenge for exploring the 
high magnetic Reynolds number Rm, the governing pa- 
rameter of the problem, because the magnetic Prandtl 
number Pm = v /rj is small, typically 10 -6 for sodium (y 
and r\ are the viscosity and magnetic resistivity). A dy- 
namo has been obtained recently within a turbulent flow 
but the high R M regime, as is the case for astrophysi- 
cal flows, remain unattainable in the laboratory. On the 



other hand, in situ observations of the Earth environ- 
ment have grown in importance recently, e.g. with the 
multi-spacecraft mission CLUSTER 3, 4]. Observations 
are quite complex but indicate clearly several features, 
such as power law energy spectra Q and intermittency 
0] (see e.g. [7j for review). One of the issues is to assess 
what kind of scaling laws obtains for both the velocity 
and the magnetic field; moreover, the flow may develop 
an anisotropic weak turbulence spectrum at small scale 
as shown recently in direct numerical simulations (DNS) 
Q and as observed in the magnetosphere of Jupiter . 

Indeed, DNS may help but remain challenging in three 
space dimensions (3D). A plethora of results concern- 
ing energy spectra in MHD have emerged recently, with 
different power laws in different regions of parameter 
space, although the boundaries between these regions 
are not fully understood and more exploration remains 
to be done. However, whatever the inertial index of the 
spectrum, one may ask whether, for correlation functions 
of higher order, similarities between hydrodynamic and 
MHD turbulence persist. It is already known that it does 
not in two space dimensions (2D) [T(| [ll], |T3] , MHD be- 
ing more intermittent than neutral fluids but the data in 
3D remains scarce. 

A further problem concerns the dissipation of energy in 
the limit of high Reynolds number R e . Mathematically, 
this is an open problem in 3D for fluids and MHD, and 
yet it is central for astrophysics where dissipative struc- 
tures, reconnection and acceleration of particles are well 
observed 0]. Intermittency (as measured by anomalous 
exponents of structure functions) and singular behavior 
are linked since the latter (except for a thin boundary 
layer delimiting the thickness of the structure) is likely to 
occur on a set of strong small-scale fluid elements highly 
localized spatially, be it vortex filaments or current and 
vorticity sheets. We thus propose in this paper an assess- 
ment of dissipation, small scale structures, intermittency, 
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TABLE I: Runs, linear resolution N, viscosity v and mag- 
netic diffusivity 77, Reynolds number Re, and Taylor Reynolds 
number R\ at peak of dissipation. 



Run 


N 


v — rj 


Re 


Rx 


I 


64 


8 x l(T ;i 


390 


180 


II 


128 


3 x 10" 3 


790 


280 


III 


256 


1.25 x 1(T 3 


1600 


430 


IV 


512 


6 x 10" 4 


3100 


630 


V 


1536 


2 x 10" 4 


10500 


1180 



and scaling laws, by analyzing a flow computed up to a 
grid resolution of 1536 3 points. 

The incompressible MHD equations read: 
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with v the velocity, b the magnetic field, V the pressure, 
pa = 1 the (uniform) density, j = V x b the current 
density, and V • v = V • b = 0. When v = r) = 0, the 
energy E = (v 2 + fc 2 ) /2, magnetic helicity Hb = (A • b) 
(with A the vector potential such that b = V x A), and 
cross helicity He = (v • b) /2, are conserved. We solve 
Eqs. J]} and ([2|) in a 3D box using periodic boundary 
conditions and a pseudospectral method dealiased by the 
standard 2/3 rule; k m - m — 1 for a box of length Lq — 2ir, 
and N regularly spaced grid points lead to a maximum 
wavenumber fc max = A^/3. At all times, we preserve 
fcu/femax < 1, where fco is the dissipation wavenumber. 

The initial conditions are constructed from a superpo- 
sition of Beltrami flows from wavenumbers k = 1 to 3, 
to which smaller-scale random fluctuations with a spec- 
trum ~ fc -3 exp[— 2(fc/fc )] 2 for k > 3 are added (see 
fj~3|). The phases of the modes with k > 3 are cho- 
sen from a Gaussian random number generator in such 
a way that the initial cross-correlation of the two fields 
is negligible: initially, Ey = Em = 0.5, He ~ 10~ 4 , 
and Hm ~ 0.45. Resolutions of runs described in this 
paper range from ./V = 64 to N = 1536 (see Table fl}. 
The largest resolution run is stopped close to the peak 
of dissipation, t = 3.7; its initial quasi-ideal phase is de- 
scribed in [l3| and the total energy spectra that develop, 
together with the ensuing anisotropy of the small scales, 
is given in @. Near the peak of dissipation, the Reynolds 
number based on the integral scale of the flow velocity is 
R e = UL/v k, 9200, and that based on the Taylor scale 
is R\ = U\/v w 1700; U is the r.m.s. velocity, the inte- 
gral scale is defined as L — 2irE~ 1 J k~ 1 E(k)dk and the 
Taylor scale as A = 2n(E/ J k 2 E(k)dk) 1 / 2 , with E(k) the 
total energy spectra. 

We focus on the fully developed turbulent regime close 
to the peak of dissipation. Figure [T] gives the varia- 
tion of the maximum of the total energy dissipation rate 
e = v (oj 2 ) + 77 (j 2 ) with R\ (w = V x v is the vor- 
ticity) for the runs of Table HI For large R\ , e seems to 
become independent of R\. This result is not entirely un- 
expected. On the one hand, the dissipation of energy is 




FIG. 1: Left: Total energy dissipation rate e as a function of 
Taylor Reynolds number R\ for several runs with the same 
initial conditions. For large R\, e seems to be independent of 
R\. Right: R\ as a function of R e calculated at the peak of 

dissipation for the same runs. The straight line indicates the 

1/2 

classical turbulent scaling R\ ~ R e . 



known to tend to a constant in the case of neutral fluids 
(b = 0) and when restricting the MHD dynamics 
to 2D (which, to lowest order, is the evolution that is 
expected in the presence of a strong uniform magnetic 
field), the e nerg y dissipation was shown similarly to be 
constant [Icllllj ; in 3D MHD, an indication that this may 
be the case as well was obtained for the Orszag-Tang vor- 
tex [l3[ although at a lower resolution (512 3 points) and 
Reynolds number (i? e ~ 5600) which did not allow for a 
clear scaling. Here, it appears that we have reached the 
beginning of an asymptotic regime where dissipation is 
constant and the Taylor Reynolds number scales as the 
square root of the Reynolds number R e (see Fig. [1]) , as 
expected for a fully developed turbulent flow. 

For e to remain constant with vanishing viscosity and 
resistivity, one can think of several scenarios; either we 
have intense dissipative structures that are more space- 
filling as R e grows, or else the structures remain sparse 
but become very sharp. Both may be happening, with 
a myriad of current sheets of intermediate to large in- 
tensity, and a few very sharp structures. When plotting 
the histogram of one component of the current intensity 
(not shown), one observes that, as the Reynolds number 
increases, the wings of the PDF stabilize at intermedi- 
ate values but substantially higher extrema are reached. 
Fig. [5] gives a 3D rendering of the current density in a 
slice of the entire domain, and in a subvolume showing 
folding and rolling of the current sheet. Visualizations of 
the time evolution of these structures confirm that the 
rolling takes place as the result of a Kelvin-Hclmholtz- 
like instability as observed in the solar wind Q ■ 

One way to determine the statistics of such structures 
is to examine the behavior of structure functions; at or- 
der p for a field u, they are defined as (|<5it||(Z)| p ), with 
6u\\ (I) = U|| (x + 1) — u || (x) with homogeneity and isotropy 
assumed and with tty the longitudinal component of the 
vector u that projects along 1. Assuming self similarity 
leads to (|<5m||(^)| p ) ~ K", with Q = ap for a scale in- 
variant (non-intermittent) field (a = 1/3 for Kolmogorov 
scaling, a = 1/4 for Iroshnikov-Kraichnan scaling). De- 
partures from such a linear scaling arc observed experi- 
mentally, obscrvationally and numerically but a normal 
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FIG. 2: (Color online) Current density in a slice of the full box 
(left), and in a subregion (right) showing folding and rolling 
of the current sheets. Vorticity organizes in the same fashion, 
although current sheets are thinner. The high intensity (dark, 
green and magenta) is concentrated in a thin layer within 
subvolumes of the flow. 

(linear) scaling occurs for third-order functions, express- 
ing the conservation laws of the ideal case: total energy 
and cross-correlation [l5|], as well as magnetic helicity 
[T^ |. In terms of the Elsasser variables = v ± b, the 
first two conservation laws lead to 

(5zJ(l)\5z±(l)\ 2 )=~e ± l, (3) 

where are the dissipation rates of . From these 
expressions, the flux of total energy e and of the cross 
correlation between the fields ec can be computed as 
a function of the scale I. The relations given by Eq. 
© as evaluated directly from the 1536 3 data near the 
peak of dissipation are shown in Fig. [3] A linear de- 
pendence with I is observed in a range of scales for both 
flux functions although the scaling is slightly better for 
the e~ flux; as a result, this is the quantity we will use 
for the extended self-similarity (ESS) analysis [ij}: in 
fluid turbulence, it is a common practice to plot struc- 
ture functions in terms of each other, the third order one 
being particularly relevant since it is proportional to I and 
can be used to define the inertial range and to improve 
the estimation of the scaling exponents. We thus deter- 
mined the anomalous scaling exponents for MHD for the 
Elsasser variables [1] using ESS. We show here the de- 
termination of these exponents for the velocity and mag- 
netic field, assuming isotropy as before. A measurable 
difference is obtained, as observed recently in the solar 
wind [To ]: it corresponds to a steeper magnetic energy 
spectrum (close to Kolmogorov scaling) and a shallower 
kinetic energy spectrum (close to Iroshnikov-Kraichnan 
scaling). Indeed, for the second order scaling exponent 
of the velocity field (v = 0-55 ± 0.01, and for the mag- 
netic field = 0.64 ± 0.01. These exponents in turn 
lead to a kinetic energy spectrum E v (k) ~ fc -1 - 55 and a 
magnetic energy spectrum Emik) ~ fc -164 . Note that 
for both fields, ( 3 ^ 1 indicating already at third order 
a departure from Kolmogorov phenomenology, and that 
the intcrmittency is stronger. However, in this simulation 
the exponents of the Elsasser variables z^ are closer to 
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FIG. 3: Top: - (^6zf(l)\8z^(l)\ 2 j as a function of the dis- 
placement / in the fully developed turbulent regime of the 
1536 3 run. A slope of 1 is indicated as a reference. Bottom: 
Scaling exponents £ p for the velocity and magnetic fields. The 
Iroshnikov-Kraichnan (IK) and Kolmogorov predictions are 
shown as a reference. Although is closer to Kolmogorov 
scaling, is far from the hydrodynamic value of 1. Note the 
measurable difference between both sets of exponents. 

Iroshnikov-Kraichnan scaling than to Kolmogorov (with 
the second order exponent near 0.6 Q because of in- 
termittency corrections). Note that the different expo- 
nents that have been observed in the solar wind [18[ 
are on the average E(k) ~ fc~ 1,6 , E v (k) ~ k~ , and 
Eb(k) ~ A: -166 [18( 1 : in solar active regions, variations 
have also been measured with a monotone decrease of 
the exponent at a given order when the strength of the 
flare augments (from Ml to XI) (l9j . 

The different scaling of the velocity and the magnetic 
field can thus be explained in terms of the different intcr- 
mittency properties of each field. Indeed, in MHD turbu- 
lence current sheets are thinner than vortex structures, a 
property that results in faster dissipation of magnetic en- 
ergy than of mechanical energy. The development of thin 
structures in the current in turn leads to a steeper spec- 
trum for the magnetic field than for the velocity field. 
Other scaling laws arise in the flow, specially at very 
high Reynolds number (run V), that in some cases have 
been previously reported in observations or predicted us- 
ing theoretical arguments. Figure 0] shows the residual 
energy spectrum En(k) = E^(k) — E v (k) [2(j, with for 
a Iroshnikov-Kraichnan scaling for total energy should 
scale as k~ 2 and for Kolmogorov scaling goes as fc -7 / 3 . 
The residual energy spectrum in the simulation is consis- 
tent with k~ 2 . The magnetic helicity Ht seems to follow 
a fc~ 10 / 3 spectrum. This scaling has also been observed in 
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FIG. 4: Residual energy En(k), magnetic helicity Ht(k) and 
ratio R(k) (see text) at the peak of dissipation in run V. 



the inverse cascade range of the magnetic helicity |2l| and 
is not well understood. It could result from the Alfvenic 
balance between Eb(k)/E v (k) and k 2 Hb/H v [2l[ where 
Hv = (v • <*>) is the kinetic helicity and where the factor 
k 2 follows from dimensional reasons; Fig. [J] also shows 

R(k) = [£ 6 (fc)/E„(fc)p 2 ff 6 (AO/#t;(fc)] -1 and although 
such a balance is plausible, a slow increase of R(k) with 
k can also be observed. 

The measurement of the energy input e (and ensu- 
ing heating) in magnctosphcric plasmas is an outstand- 
ing problem that the present CLUSTER mission helps 
unravel and that future missions, such as the Magnc- 
tospheric MultiScale (MMS) to be launched in 2013 is 
designed to study. From solar wind observations, it can 
be measured by using the exact scaling laws used here, 
which seem to be reasonably followed by the numerical 
data. The fact that the present study shows the con- 
stancy of e with Reynolds number indicates that energy 
is transferred efficiently to small scales in MHD as long 
as sufficient scale separation is available. Kinetic plasmas 



effects will come into play as the cascade meets, e.g., the 
ion-cyclotron frequency, leaving open the issue of what 
follows at smaller scale, but the present results imply 
that energy can be cascaded rapidly (independently of 
the Reynolds number) to the smallest available scales by 
MHD turbulence. Moreover, dissipation is achieved in 
localized regions with strong magnetic field gradients, in 
the form of current sheets. These extreme events, more 
probable at small scales than what is expected from a 
normal distribution, represent a break down of scale in- 
variance and give rise to intermittency. The thin cur- 
rent sheets result in a more intermittent magnetic field 
than velocity field, and in turn make the magnetic energy 
spectrum steeper than the kinetic energy spectrum. Re- 
markably, the second order scaling exponents and spec- 
tral indices for the kinetic, magnetic, and total energy in 
the MHD simulation at largest Reynolds number are in 
good agreement with the ones reported for the solar wind. 
However, care must be taken when extracting conclusions 
about scaling laws in MHD turbulence. Simulations have 
been reported where the total energy spectrum follows 
different power laws depending on properties of the forc- 
ing (22| , and in the solar wind variations in the total 
energy spectrum from ~ k~ 3 / 2 to ~ k~ 5 ^ 3 have been ob- 
served [l8[ . It is not our intention to say that MHD tur- 
bulence has unique scaling properties represented by our 
simulations, but rather that the determination of scal- 
ing laws in MHD turbulence, and the explanation of the 
results from solar wind observations, require the study 
of often neglected phenomena as intermittency, and the 
measurement of high order statistics of the velocity and 
the magnetic fields. 

The largest run was performed through a BTS grant at 
NCAR, sponsored by NSF. PDM acknowledges support 
from grant UBACYT X4-68/08 and from the Carrera del 
Investigador Cientifico of CONICET. Three-dimensional 
visualizations use VAPOR, a software for interactive vi- 
sualization and analysis of terascale datasets W^J . 
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